1 



On feasibility, stability and performance in 
distributed model predictive control 

Pontus Giselsson and Anders Rantzer 

Department of Automatic Control 

Lund University 
Box 118, SE-221 00 Lund, Sweden 

{pontusg,rantzer}@oontrol.lth.se 



o 

CD 
00 



U 
O 



> 

O 

m 



X 



Abstract — In distributed model predictive control (DMPC), 
where a centralized optimization problem is solved in distributed 
fashion using dual decomposition, it is important to keep the 
number of iterations in the solution algorithm, i.e. the amount 
of communication between subsystems, as small as possible. 
At the same time, the number of iterations must be enough 
to give a feasible solution to the optimization problem and to 
guarantee stability of the closed loop system. In this paper, 
a stopping condition to the distributed optimization algorithm 
that guarantees these properties, is presented. The stopping 
condition is based on two theoretical contributions. First, since 
the optimization problem is solved using dual decomposition, 
standard techniques to prove stability in model predictive control 
(MFC), i.e. with a terminal cost and a terminal constraint set 
that involve all state variables, do not apply. For the case without 
a terminal cost or a terminal constraint set, we present a new 
method to quantify the control horizon needed to ensure stability 
and a prespecified performance. Second, the stopping condition 
is based on a novel adaptive constraint tightening approach. 
Using this adaptive constraint tightening approach, we guarantee 
that a primal feasible solution to the optimization problem is 
found and that closed loop stability and performance is obtained. 
Numerical examples show that the number of iterations needed 
to guarantee feasibility of the optimization problem, stability 
and a prespecified performance of the closed-loop system can 
be reduced significantly using the proposed stopping condition. 

Index Terms — Distributed model predictive control, perfor- 
mance guarantee, stability, feasibility 



I. Introduction 

Distributed model predictive control (DMPC) can be divided 
into two main categories. In the first category, local opti- 
mization problems that are solved sequentially and that take 
neighboring interaction and solutions into account, are solved 
in each subsystem. This is done in 1 1 1 for linear systems and in 
ij^l for nonlinear systems. In Q a DMPC scheme is presented 
in which stability is proven by adding a constraint to the 
optimization problem that requires a reduction of an explicit 
control Lyapunov function. In ||4|, 15] stability is guaranteed 
for systems satisfying a certain matching condition and if the 
coupling interaction is small enough. In the second category, 
to which the current paper belong, a centralized optimization 
problem with a sparse structure is solved using a distributed 
optimization algorithm. This approach is taken in ||6] where 
stability is guaranteed in every algorithm iteration. A drawback 
to this method is that full model knowledge is assumed in 
each node. Other approaches in the DMPC literature rely on 
dual decomposition to solve the centralized MPC problem in 



distributed fashion. This approach is taken in, e.g. ff], fSl, 
|9|, where a (sub)gradient algorithm is used to solve the dual 
problem and in [iTOl where the algorithm is based on the 
smoothing technique presented in fTl |. Among these, the only 
stability proof is given in [12|, |9|, where a terminal point 
constraint is set to the origin, which is very restrictive. 

One reason for the lack of stability results in DMPC based 
on dual decomposition, is that the standard techniques to prove 
stability in MPC do not apply. In MPC, terminal costs and 
terminal constraint sets that involve all state variables are used 
to show stability of the closed loop system, see |[T3l , lfT4ll . 
This is not compatible with dual decomposition. However, 
results for stability in MPC without a terminal constraint 
set or a terminal set, which fits also the DMPC framework 
used here, are available ||T5l . |fT6l . In ||T6I . a method to 
quantify the minimal control horizon that guarantees stability 
and a prespecified performance is presented. This is based on 
relaxed dynamic programming IfTTI . llTSi and a controllability 
assumption on the stage costs. In the current paper, we take 
a similar approach to quantify the control horizon needed 
to guarantee stability and a prespecified performance. The 
advantages of our approach over the one in [16] are twofold; 
we can, by solving a mixed integer linear program (MILP), 
verify our controllability assumption, further we get an explicit 
expression that relates the parameter in the controllability 
assumption with the obtained closed loop performance. 

Besides the stability result, the main contribution of this 
paper is a stopping condition for DMPC controllers that use 
a distributed optimization algorithm based on dual decompo- 
sition. We use the distributed algorithm presented in |19 |, but 
any duality-based distributed algorithm, such as the standard 
dual ascent or ADMM |20|, can be used. These duality 
based algorithms suffer from that primal feasibility is only 
guaranteed in the limit of iterations. Constraint tightening, 
which was originally proposed for robust MPC in [21], can 
also be used to generate feasible solutions within finite number 
of iterations, see |22|. However, the introduction of constraint 
tightening complicates stability analysis since the optimal 
value function without constraint tightening is used to show 
stability, while the optimization is performed with tightened 
constraints. This problem is addressed in |22 | by assuming 
that the difference between the optimal value functions with 
and without constraint tightening is bounded by a constant. 
However, to actually compute such a constant is very difficult. 
The stopping condition in this paper is based on a novel adap- 
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tive constraint tightening approach that ensures feasibility w.rt. 
the original constraint set with a finite number of algorithm 
iterations. In addition, the amount of constraint tightening 
is adapted until the difference between the optimal value 
functions with and without constraint tightening is bounded 
by a certain amount. This adaptation makes it possible to 
guarantee, besides feasibility of the optimization problem, also 
stability of the closed-loop system, without stating additional, 
unquantifiable assumptions. 

The paper is organized as follows. In Section |ll]we introduce 
the problem and present the distributed optimization algorithm 
in |19|. In Section Hill the stopping condition is presented and 
feasibility, stability, and performance is analyzed. Section |IV] 
is devoted to computation of a controllability parameter in the 
controllability assumption. A numerical example that shows 
the efficiency of the proposed stopping condition, is presented 
in Section |V] Finally, in Section [VT] we conclude the paper 

II. Problem setup and preliminaries 
We consider linear dynamical systems of the form 

xt+i = Axt + But, xq ^ X (1) 

where Xt € M" and ut G M™ denote the state and control 
vectors at time t and the pair {A, B) is assumed controllable. 
We introduce the following state and control variable partitions 

xt^[{x]rMf,...,{xrrr, 



where the local variables x\ e 



and u\ e 



(2) 
(3) 
The A and 



B matrices are partitioned accordingly 

/ All ■ ■ ■ Aim \ ( Bu 
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MM / 
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These matrices are assumed to have a sparse structure, i.e., 
some Aij — and Bij — and the neighboring interaction is 
defined by the following sets 

M = {j e {1, . . . , M} I if Ay ^ or By ^ 0} 

for i — 1, . . . , M. This gives the following local dynamics 



E 



A. 



Bijul 



for i = 1, . . . , A/. The local control and state variables are 
constrained, i.e., G Ui and a;' S Xi. The constraint sets, Xi, 
Ui are assumed to be bounded polytopes containing zero in 
their respective interiors and can hence be represented as 



X, = {x' 



where CI G 



I cw < <} 



and dl G 



R^q"' . We also denote the total number of linear inequalities 
describing all constraint sets by ric := X^tf^i ("c i + ric j). 
The global constraint sets are defined from the local ones 
through 



X =XiX ...X X, 



M, 



U =Ui X . . . X Um- 



We use a separable quadratic stage cost 



u) 



M /M \ 

Y,h{x\u') = - Y.{xYQ,x' + {u'fRy 

1=1 \i=l / 

where Qi e §"!^_ and Ri e for z = 1,...,A/ and 

§"_,_ denotes the set of symmetric positive definite matrices 
in M"^". The optimal infinite horizon cost from initial state 
a; G A" is defined by 



Vo^{x):= min yj{xuUt) (4) 

s.t. Xt ^ X , Ut eU 

xt+i = Axt + But 

Xq = X. 

Such infinite horizon optimization problems are in general 
intractable to solve exactly. A common approach is to solve 
the problem approximately in receding horizon fashion. To this 
end, we introduce the predicted state and control sequences 
{Zr}^SQ and {vT}r=o '■^^ corresponding stacked vectors 



[Zq , . . . , zjf_i]'^ , 



r T T iT 



where Zr and are predicted states and controls t time steps 
ahead. The predicted state and control variables Zr, Vr are 
partitioned into local variables as in dU and (|3]l respectively. 
We also introduce the following stacked local vectors 



[K)^■••,(^]v-l) 



TiT 



[K)^...,(^;]v-l) 



TiT 



Further, we introduce the tightened state and control constraint 
sets 

(1 - &)X, = {x' G K"' I Clx' < (1 - 5)4}, 
(1 - 5)U^ = {u' G M"' I < (1 - 

where 5 G (0,1) decides the amount of relative constraint 
tightening. The following optimization problem, which has 
neither a terminal cost nor a terminal constraint set, is solved 
in the DMPC controller for the current state x G M" 

Af-l 



Vt{x) := min i(zr,Vr) 



T = 



S.t. z^e{l-S)X,T^O,...,N-l 
Vr G {l-S)U, r = 0, ...,iV-l 
Zr+i = AZr + BVr, T ^Q,...,N -2 
zo = X. 



By stacking all decision variables into one vector 



(5) 



(6) 



y — [Zq , . . . , Z]s^_i,Vq , 

the optimization problem (|5]l can more compactly be written 
as 



V^{x):= min iy^Hy 
y 

s.t. Ay = hx 

Cy < (1 - S)d 



(7) 



where H G §^'+"''^,A G M«(w-i)x(n+m)JV^ ^ ^ 
accordingly. Such sparse optimization problems can be solved 



^>0 



are built 
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in distributed fashion using, e.g., the classical dual ascent, the 
alternating direction of multipliers method (ADMM) [20], or 
the recently developed algorithm in 119|. The algorithm in 
|fT9l is a dual accelerated gradient algorithm and is used in 
the current paper for simplicity. Distribution of these methods 
are enabled by solving the dual problem to (|5]l. 

The dual problem to Q is created by introducing dual 
variables A G for the equality constraints and dual 

variables /i, e ]R>q " for the inequality constraints. As shown 
in ||T9l , the dual problem can explicitly be written as 

max --{A^X + C'^ nfli-^{A^ X + C^/x)- 

A,/i.>0 2 

- X^hx- fjFd{l- S) (8) 

and we define the minimand in ^ as the dual function for 
initial condition x £ M", i.e., 

D%ix,X,n) := -i(A^A + C^/x)^H-i(A^A + C^/x)- 

~ X^hx- fi^d{l~ 6). (9) 

The distributed algorithm presented in [19] that solves (|7), is 
a dual accelerated gradient method described by the following 
global iterations 

y'' = -H"\A^A'= + C^/x'') (10) 

K + 2 

^^k^ '^{X" - X'-') + yiAy'^ - hx) (12) 
K + 2 L 

+ l(Cy^-d(l-<5))) (13) 

where k is the iteration number and L = 
II [A^, C^]'^H^i[A^, C"^]]], which is the Lipschitz constant 
to the gradient of the dual function (|9]l. The reader is referred 
to IIT9I for details on how to distribute the algorithm (fT0i- (fT3] l. 

A. Notation 

We define by N>t the set of natural numbers t > T. 
The norm || • || refers to the Euclidean norm or the in- 
duced Euclidean norm unless otherwise is specified and (•, •) 
refers to the inner product in Euclidean space. The norm 
||x||t = V x^Tx. The optimal state and control sequences 
to Q for initial value x and constraint tightening 5 are 
denoted {z*{x,5)}'^~q and {v*{x,5)}^^q respectively and 
the optimal solution to the equivalent problem y*(x,5). 
The state and control sequences for iteration k in (fT0]l-(fT3li 
are denoted {z^ix^SYl^ZQ and {w^(a;, (5)}^Jg^ respectively. 
The initial state and constraint tightening arguments (a;, 5) are 
dropped when no ambiguities can arise. 

B. Definitions and assumptions 

We adopt the convention that V^{x) — 00 for states x e K" 
that result in d?) being infeasible. We define by Xqo the set 



for which (|4]i is feasible and we define the minimum of the 
stage-cost £ for fixed x as 

£*(x) :=min^(a;, 7i) = —x'^Qx. 

u<£U 2 

Further, k is the smallest scalar such that kQ — A^QA > 0. 
The state sequence resulting from applying {wr}^=o^ to ([U is 
denoted by {5,t}t=o^ l-^-' 

£,r+l= A£,r+BVr, io=X. (14) 

We introduce ^ — [{io)^ , ■ ■ ■ t {^n -li^Y' ^^'^ define the 
primal cost 

e{£.r, Vr) if I € A-^ and v 

and ([141l holds 
oo else 

(15) 

where and are the state and control constraints for the 
full horizon. We also introduce the shifted control sequence 
V, = [iv,)^,...,ivN-iV,0'^f. We have Pjv(x,v^) > 
Va,(x) and Pn {Ax + Bv^,v'^) > Vn{Ax + Bv^) for every al- 
gorithm iteration k. We denote by the state sequence 
that satisfies (fl4] i using controls {w^'}^q^. The definition of 
the cost ( fTsT i implies 

Pjv(x,v'=) = Pn{Ax + Bv^,v':) + £{x,v^o) - £*iMN-i) 

(16) 

if eU, xeX and A^%_^ e X. 

III. Stopping condition 

Rather than finding the optimal solution in each time step 
in the MFC controller, the most important task is to find 
a control action that gives desirable closed loop properties 
such as stability, feasibility, and a desired performance. Such 
properties can sometimes be ensured well before convergence 
to the optimal solution. To benefit from this observation, a 
stopping condition is developed that allows the iterations to 
stop when the desired performance, stability, and feasibility 
can be guaranteed. Before the stopping condition is introduced, 
we briefly go through the main ideas below. 

A. Main ideas 

The distributed nature of the optimization algorithm makes 
it unsuitable for centralized terminal costs and terminal con- 
straints. Thus, stability and performance need to be ensured 
without these constructions. We define the following infinite 
horizon performance for feedback control law i/ 

oo 
t=0 

where Xf+i — Axt + Bi'{xt) and xq — x. For a given 
performance parameter a € (0, 1] and control law i/, it is 
known (cf. ITtI . ifTSl . lfT6l ) that the following decrease in the 
optimal value function 

Vj^ixt) > V^iAxt + Bv{xt)) + al{x, v{xt)) (18) 
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for every t G N>o gives stability and closed loop performance 
according to 

c^V^i^) < Voo{x). (19) 

Analysis of the control horizon N needed for an MPC con- 
trol law without terminal cost and terminal constraints such 
that ( fTsT i holds, is performed in fTSl, fTSl and also in this 
paper. Once a control horizon N is known such that ( fTSl l is 
guaranteed, the performance result (|T9] relies on computation 
of the optimal solution to the MPC optimization problem 
in every time step. An exact optimal solution cannot be 
computed and the idea behind this paper is to develop stopping 
conditions that enable early termination of the optimization 
algorithm with maintained feasibility, stability, and perfor- 
mance guarantees. The idea behind our stopping condition 
is to compute a lower bound to V^{x) through the dual 
function D^j^{x, A'", ^jl^) and an upper bound to the next step 
value function V^{Ax + Bv^) through a feasible solution 
Pn{Ax + i?WQ,v^). If at iteration k the following test is 
satisfied 

D%{x, A^ /x'^) > Pjv(At + Bv^^y^) + al{x, «o^) (20) 

the performance condition ( fTSl l holds since 

VU^) > D%(x, A^ m'^) > PNiAx + v^-) + aiix, v^) 
> V^{Ax + Bv^)+ae{x,v^). 

This implies that stability and the performance result (|T9^ 
can be guaranteed with finite algorithm iterations k by using 
control action Vq. 

The test (|20] | includes computation of Pn{Ax + Bvqjv'^) 
which is a feasible solution to the optimization problem in the 
following step. A feasible solution cannot be expected with 
finite number of iterations k for duality-based methods since 
primal feasibility is only guaranteed in the limit of iterations. 
Therefore we introduce tightened state and control constraint 
sets (1 - S)X, (1 - 5)U with d £ (0,1) and use these in 
the optimization problem. By generating a state trajectory 
{S,t}t=o from the control trajectory {v'^}!^Sq that satisfies the 
equality constraints (fl4l l. we will see that {S,'^}!^Sq satisfies the 
original inequality constraints with finite number of iterations. 
Thus, a primal feasible solution Pn{Ax + Bvq,v'^) can be 
generated after a finite number of algorithm iterations k. 
However, since the optimization now is performed over a 
tightened constraint set, the dual function value _D^(x, A, /x) 
is not a lower bound to V^{x) and cannot be used directly in 
the test ( |20] | to ensure stability and the performance specified 
by (T% . In the following lemma we show a relation between 
the dual function value when using the tightened constraint 
sets and the optimal value function when using the original 
constraint sets. 

Lemma 1: For every x £ K", A e M"(^-i) and ^l e M>^<^ 
we have that 

V^{x) > D%{x,\ti) - Sfi^d. 

Proof. From the definition of the dual function (|9|l we get that 
D%{x, A, /x) = D%{x, A, fi) + Sd'^fi. 



By weak duality we get 

V^ix) > D%{x, A, /x) = D%ix, A, fi) - M^/x. 

This completes the proof. □ 

The presented lemma enables computation of a lower bound 
to V^{x) at algorithm iteration k that depends on Sfi^d. By 
adapting the amount of constraint tightening 6 to satisfy 

SifJ.'^fd < e£*{x) (21) 

for some e > and use this together with the following test 

D%{x,\^,n'') > Pn{Ax + Bv^,^^) + at{x,v^) (22) 

we get from Lemma [T] and if ( 1211 1 and ( l22l l holds that 

v^{x)>D%{x,\\^,'^)~5{^l^fd 

> Pn{Ax + Bv'^, v^) + al{x, v^) - er{x) 

> V^iAx + Bv^) + {a- e)i{x,v^). 

This is condition ( fTSl l, which guarantees stability and perfor- 
mance specified by (fT9] l if a > e. 

B. The stopping condition 

Below we state the stopping condition, whereafter parameter 
settings are discussed. 

Algorithm 1: Stopping condition 

Input: X 

Set: k = 0,1 = 0,5 = 5n,it 

Initialize algorithm (fTOli-lflJTl with: 

A° = A"^ = 0, /xO = = and y" = y-^ = 0. 

Do 

\tD%{x,\\tJ.^)>PN{x,^^)-j^t{x) 
or Sd^n*' > ei*{x) 

Set (5 ^ (5/2 // reduce constraint tightening 

Set ; ^ / + 1 

Set fc = // reset step size and iteration counter 

End 

Run Ak iterations of (fT0]l-(fT3]l 

Set fc ^ fc + AA; 
Until D%{x,\^,n^) > Pn{Ax + Bv^,v';) + a£{x,v^) and 

Sd'^fi'' < e£*{x) 
Output: Vq 



In Algorithm [T] four parameters need to be set. The first is 
the performance parameter a G (0, 1] which guarantees closed 
loop performance as specified by ([19). The larger a, the better 
performance is guaranteed but a longer control horizon N 
will be needed to guarantee the specified performance. The 
second parameter is an initial constraint tightening parameter, 
which we denote by (5i„it G (0, 1], from which the constraint 
tightening parameter 6 will be adapted (reduced), to satisfy 
([2Tl l. A generic value that always works is S-mit = 0.2, i.e., 
20% initial constraint tightening. The third parameter is the 
relative optimality tolerance e > where e < a. The e must 
be chosen to satisfy (1251 ). Finally, Afc, which is the number 
of algorithm iterations between every stopping condition test. 
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should be set to a positive integer, typically in the range 5 to 
20. 

Except for the initial condition x. Algorithm [1] is always 
identically initialized and follows a deterministic scheme. 
Thus, for fixed initial condition the same control action is 
always computed. This implies that Algorithm [T] defines a 
static feedback control law, which we denote by j/jv- We get 
the following closed loop dynamics 



xt+i = Axt + Bi/Nixt), 



Xo 



The objective of this section is to present a theorem stating 
that the feedback control law function i^n satisfies dom(i/iv) 2 
int(X^), where 



{x G R" I V^{x) < oo and Az*n_^{x,0) G int(A')} 

(23) 



C X'^ for Si > 62- First, however we state 



which satisfies 
the following definition. 

Definition 1: The constant $Ar is the smallest constant such 
that the optimal solution {z*{x, 0)}f^o^ {^r (^^ ^)}t=o ^° © 
for every x G X^ satisfies 

r(z^_i(5,0)) < $jv^(x,«o*(S,0)) 

for the chosen control horizon N . 

The parameter $jv is a measure that compares the first and 
last stage costs in the horizon. In Section |IV] a method to 
compute $jv is presented. 

Remark 1: In ifTsl . lfT6l an exponential controllability on 
the stage costs is assumed, i.e., that for C > 1 and a G (0, 1) 
the following holds for t = 0, . . . , N — 1 

e*{z*{x,0),v*{x,0)) < Ca''e{x,v*(x,0)). 

This impHes <^>n < Ca^-^. 

We also need the following lemmas, that are proven in Ap- 
pendixTAl AppendixFBl and Appendi^jTCl respectivelv. to prove 
the upcoming theorem. 

Lemma 2: Suppose that e > and 5 G (0,1]. For every 
X G X'L we have for some finite k that 



D%{x,X\^l'')>PN{x,w^) 



ei*ix). 



(24) 



1 






'z;{x,o) 




2 






v;{x,o) 





Lemma 3: Suppose that e > and 5 G (0,1]. For every 
X G X'^ and algorithm iteration k such that ( l24l i holds we 
have for r = 0, . . . , - 1 that 

2 

<er{x)+d{n^fd 

H 

where H = blkdiag(Q, R). 

Lemma 4: Suppose that e > and 5 G (0, 1]. For x G X^J^ 
but X ^ X^ we have that 5{h^YA > e£*(x) with finite k. 

We are now ready to state the following theorem, which is 
proven in Appendi?(GD] 

Theorem 1: Assume that e > 0, (Jjnit G (0, 1] and 



a < 1 



-K(V2e+ v/$Jv)'(V2e + l)^ (25) 

Then the feedback control law v^q, defined by Algorithm [T] 
satisfies dom(i/Ar) I) int(X5^). Further 

V^{x) > V^iAx + Bvn{x)) + (a - e)€(x, vn^x)) (26) 



holds for every x G dom(i^N). 

Corollary 1: Suppose that a < 1 — k^^t and that v^{x) = 
vl{x,Q). Then 

V^{x) > V^{Ax + Bv*^{x)) + al{x, i^^(a;)). 

holds for every x G X^^. 

Proof. For every x G X^^r we have 

N-l 

V^{x) = J2 <) + ^(^^^-1,0) - 0) 

> V^{Ax + Bv*^{x)) + l{x, «o*) - 0) 

> V^[Ax + Bv*^{x)) + i{x, v*q) - k£{z*^_^, 0) 

> V^iAx + Bv*^{x)) + (1 - k^n)1{x,v*) 

where the first inequality holds since Az^ ^^ G A" by con- 
struction of X^, the second due to the definition of k and the 
third due to the definition of . □ 

Remark 2: By setting e = in Theorem [T| we get a < 
1 — K$ jv as in Corollary [T] 

C. Feasibility, stability and performance 

The following proposition shows one-step feasibility when 
using the feedback control law i^jy. 

Proposition 1: Suppose that a satisfies ( l25T l. For every xt G 
int(X5^) we have that xt+i = Axt + Bv^ixt) G X. 
Proof. From Theorem [T] we have that xt G dom(z^jv) and 
from Algorithm [U we have that PM{xt+i,'v'l) < 00 which, by 
definition, implies that xt+i EX. □ 

The proposition shows that Xt+i is feasible if xt G dom(i^Ar). 
We define the recursively feasible set as the maximal set such 
that 

Xrf = {xeX \ Ax + Bvn{x) G Xrf} 

In the following theorem we show that Xif is the region of 
attraction and that the control law vm achieves a prespecified 
performance as specified by (fTTI i. 

Theorem 2: Suppose that a > e satisfies ( l25T l. Then for 
every initial condition x G Xif we have that ||a;f|| — > as 
t — > 00 and that the closed loop performance satisfies 

(a-e)K^~(x) < Foo(x). (27) 

Further, Xif is the region of attraction. 

Proof. From the definition of X^f we know that a; = xo G X^f 
implies xt G X^f for all t G N>o. Since, by construction, 
Xrf ^ int(X5^) C dom(i^jv) we have from Theorem [T] that 
(|2^ holds for all xt, t G N>o. In HS). Proposition 2.2] it 
was shown, using telescope summation, that (l26T l implies dZTl i. 
Further, since the stage cost £ satisfies [16, Assumption 5.1] 
we get from |16, Theorem 5.2] that ||xt|| — as t 00. 

What is left to show is that Xrf is the region of attraction. 
Denote by Xioa the region of attraction using v^. We have 
above shown that X^f C Xioa. We next show that Xroa C X^f 
by a contradiction argument to conclude that X^f ~ Xioa. 
Assume that there exist x G Xroa such that x ^ Xrf. If 
X G Xroa the closed loop state sequence {xtj^p is feasible 
in every step (and converges to the origin) and consequently 
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{Axt + Bi'N{xt)}'^Q is feasible in every step. This is exactly 
the requirement to have x G Xj-f, which is a contradiction. 
Thus Xif C Xioa ^ Xrf which implies that X^f — Xroa- 
This completes the proof. □ 

To guarantee a priori that the control law z^jv achieves 
the performance (|27i specified by a, we need to find a 
control horizon N such that the corresponding controllability 
parameter $jv satisfies (IZST l. This requires the computation of 
controllability parameter which is the topic of the next 
section. 

IV. Offline controllability verification 

The stability and performance results in Theorem |2] rely 
on Definition [T] For the results to be practically meaningful 
it must be possible to compute <i>jv in Definition [T] In this 
section we will show that this can be done by solving a 
mixed integer linear program (MILP). For desired performance 
specified by a, we get a requirement on the controllability 
parameter through (|25t for Theorem [Tjand Theorem 2]to hold. 
We denote by the largest controllability parameter such 
that Theorem [T] and Theorem |2] holds for the specified a. This 
parameter is the one that gives equality in (IZSl l, i.e., satisfies 

a = 1 - e- K(\/2i+ y$^)^(%/2^+ 1)^ (28) 

for the desired performance a and optimality tolerance e. The 
parameters a and e must be chosen such that > 0. The 
objective is to find a control horizon N such that the cor- 
responding controllability parameter satisfies <I>jv < ^a- 
First we show that for long enough control horizon N there 
exist a $Ar < 

Lemma 5: Assume that a and e are chosen such that "I>q, > 
where is implicitly defined in (|28] |. Then there exists 
control horizon N and corresponding controllability parameter 

Proof. Since Xj-f is the region of attraction we have Xrf C 
Xoo. This in turn impUes that (|7]i is feasible for every control 
horizon N € N>i due to the absence of terminal constraints. 
We have 

JV-2 

> Vn-i{x)+£{z*^_^,vI^_^). 

Since the pair (A, B) is assumed controllable and since (|7]i 
has neither terminal constraints nor terminal cost we have for 
some finite M that M > Voo{x) > Vn{x) > Vn-i{x). Thus 
the sequence {Vn {x)}]^^q is a bounded monotonic increasing 
sequence which is well known to be convergent. Thus, for 
N > N where N is large enough the difference Vn{x) — 
Vn-i{x) is arbitrarily small. Especially i{z'^_^,v'^_^) = 
i*iz*j^-i) < VNix) - VN-iix) < <i>Jix,v^) since > 0. 
That is, for long enough control horizon N > N, < 
This completes the proof. □ 

The preceding Lemma shows that there exists a control 
horizon N such that < if > for the chosen 
performance a and tolerance e. The choice of performance 
parameter a gives requirements on how e can be chosen to 



give $Q > 0. Larger e requires smaller to satisfy (1281) 
which in turn requires longer control horizons N since $jv 
must satisfy $Ar < In the following section we address 
the problem of how to compute the control horizon N and 
corresponding $Ar such that the desired performance specified 
by a can be guaranteed. 



A. Exact verification of controllability parameter 

In the following proposition we introduce an optimization 
problem that tests if the controllability parameter $jv cor- 
responding to control horizon N satisfies < for 
the desired performance specified by a. Before we state the 
proposition, the following matrices are introduced 

T = blkdiag(0, . . . , 0, -g, 0, . . . , 0, -R) 
S' = blkdiag(0, ...,0,/,0,...,0) 

where Q and R are the cost matrices for states and inputs and 
$Q is the required controllability parameter for the chosen a. 
Recalling the partitioning (|6]l of y implies that 

y^Ty ^aVoRvo - ^^-iQ^Af-i - vJf_^RvN-i 
Sy = ZN-i 

Proposition 2: Assume that "I>q, > satisfies (l28T l for the 
chosen performance parameter a and optimality tolerance e. 
Further assume that the control horizon N is such that 

- mm i ($„i^Q5 + y^Ty) (29) 

s.t. X eX% 

y = argmin V^(x) 

then $Ar < $Q,. 

Proof. First we note that x — gives y = and ^aX^Qx + 
y^Ty — 0, i.e., we have that is always a feasible solution. 
Further, (|29t implies for every x £ X^^ that 

< $aX^Qa; + y^Ty = v*o) - £(z^„i, 

= <^J{x,v*)-r{z%^,) 

since = 0. This is exactly the condition in Definition [T] 

Since $jv is the smallest such constant, we have $Ar < $q, 
for the chosen control horizon N and desired performance a 
and optimality tolerance e. □ 

The optimization problem ( [29] l is a bilevel optimization 
problem with indefinite quadratic cost (see ||231 for a survey 
on bilevel optimization). Such problems are in general NP- 
hard to solve. The problem can, however, be rewritten as an 
equivalent MILP as shown in the following proposition which 
is a straightforward application of [24. Theorem 2]. 

Proposition 3: Assume that $a satisfies (|28l l for the chosen 
performance parameter a and optimality tolerance e. If the 
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control horizon N is such that the following holds 

= min _ i (d^^c^i + _^ d^^") (30) 

s.t. /3f e {0, 1} , e {0, 1} , /3P e {0, 1} 

Upper level 

Primal and dual feasibility 

C xX — dx — ^ 
< , > 

CxASy -dx-s' =0 
[_ < , A*^^ > 
Stationarity 

$,gx + iCxV^i^^ - b^A^^^2 ^ 

AA^^i = 
[ CA^^^i - ^ 
Complementarity 

/3f = 1 ^ = , /3f = ^ fif^^ = 
/3f 1 = 1 ^ sf = , = ^ /Ltf 1 = 
_ [ /3f 2 ^ 1 ^ = , /3f 2 = ^ 2 = 

Lower level 

Primal and dual feasibility 

Ay -hx = 

Cy - d - s 
[ s<0 , >0 
Stationarity 

[ Hy + A^A^ + C^^^ = 
Complementarity 

_ [ /3f = 1 ^ = , /3,f = ^ = 

where all /3, /i, A, s and x, y are decision variables, then > 

Proof. The set can equivalently be written as 

X% = {xeW I Ay* {x, 0) = bx, Cy* (a;, 0) < d, 

CxASy*ix,0)<dx,CxX<dx}. (31) 

We express the set X^^ in ( |29] ) using ( [3T| i. The equivalence 
between the optimization problems ( |30] | and (|29] l is established 
in ll24l Theorem 2]. The remaining parts of the proposition 
follow by applying Proposition |2] □ 

The transformation from ( |29t to (|30t is done by expressing 
the lower level optimization problem in ( |29t by its sufficient 
and necessary KKT conditions to get a single level indefi- 
nite quadratic program with complementarity constraints. The 
resulting indefinite quadratic program with complementarity 
constraints can in turn be cast as a MILP to get ( [30l l. 

Remark 3: Although MILP problems are NP-hard, there are 
efficient solvers available such as CPLEX and GUROBl. There 
are also solvers available for solving the bilevel optimiza- 
tion problem (|29] l directly, e.g., the function solvebilevel in 
YALMIP |25|. 

If the chosen control horizon N is not long enough for 
£ ^a, different heuristics can be used to choose a new 
longer horizon to be verified. One heuristic is to assume 
exponential controllability as in RemarklT] i.e., that there exist 



constants C > 1 and ct € (0, 1) such that 

Ca^i{x,vl)>i{z''^,v':) (32) 

for all T = 0, . . . , — 1. The C and cr-parameters should be 
determined using the optimal solution y to O for the x that 
minimized ( [30b in the previous test. Under the assumption 
that ( l32b holds as N increases, a new guess on the control 
horizon can be computed by finding the smallest such 
that Ccr^-i < 

B. Controllability parameter estimation 

The test in Proposition [3] verifies if the control horizon 
N is long enough for the controllability assumption to hold 
for the required controllability parameter Thus, an initial 
guess on the control horizon is needed. A guaranteed lower 
bound can easily be computed by solving O for a variety 
of initial conditions x and compute the worst controllability 
parameter, denoted by <&Ar, for these sample points. If the 
estimated controllability parameter > we know that 
the control horizon need to be increased for d30t to hold. If 
instead $7v £ 'I'q the control horizon N might serve as a 
good initial guess to be verified by d30] l. 

Remark 4: For large systems, (l3Qt may be too complex to 
verify the desired performance. In such cases, the heuristic 
method mentioned above can be used in conjunction with 
an adaptive horizon scheme. The adaptive scheme keeps the 
horizon fixed for all time-steps until the controllability as- 
sumption does not hold. Then, the control horizon is increased 
to satisfy the assumption and kept at the new level until the 
controllability assumption does not hold again. Eventually the 
control horizon will be large enough for $Ar < and the 
horizon need not be increased again. 

V. Numerical example 

We evaluate the efficiency of the proposed distributed feed- 
back control law vj<! by applying it to a randomly generated 
dynamical system with sparsity structure that is specified in 
[26, Supplement A.l]. The random dynamics matrix is scaled 
such that the magnitude of the largest eigenvalue is 1.1. The 
system has 3 subsystems with 5 states and 1 input each. 
All state variables are upper and lower bounded by random 
numbers in the intervals [0.5, 1.5] and [—0.15, —0.05] respec- 
tively and all input variables are upper and lower bounded by 
random numbers in the intervals [0.5,1.5] and [-0.5,-1.5] 
respectively. The stage cost is chosen to be 

^i^^i^'^i) — ~t~ Xli 

for i — 1,2,3. The suboptimality parameter is chosen 
a = 0.01. According to Theorem [T] to quantify the control 
horizon N{a), the optimality tolerance e must be chosen 
and K computed, where k is the smallest constant such that 
kQ >z A^QA. We get k = 1.22 and choose e = 0.005. Using 
( l25b . we get $Ar(o.oi) ^ 0.51. Verification by solving the 
MILP in ( [30l ) gives that the smallest control horizon N{0.01) 
that satisfies $Ar(o.oi) < 0.51 is iV(O.Ol) = 6. 

Table U presents the results. The first column specifies 
the stopping condition used, "stop, cond." for the stopping 
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TABLE I 

Experimental results for different performance requirements 
a and different initial constraint tightenings h-^-aw,- 



Algorithm comparison, ol = 0.01, N = 6 



condition 


e 


"Jinit 


avg. # iters 


max # iters 


avg. 5 


stop. cond. 


0.005 


0.001 


288.3 


506 


0.001 


stop. cond. 


0.005 


0.01 


151.5 


260 


0.01 


stop. cond. 


0.005 


0.05 


73.7 


237 


0.049 


stop. cond. 


0.005 


0.1 


70.7 


236 


0.057 


stop. cond. 


0.005 


0.2 


72.8 


236 


0.060 


stop. cond. 


0.005 


0.5 


69.2 


234 


0.076 


opt. cond. 


0.005 


0.001 


324.5 


506 


0.001 


opt. cond. 


0.005 


0.01 


171.5 


260 


0.01 



condition presented in Algorithm [T] and "opt. cond." for a 
optimality conditions. The second column specifies the duality 
gap tolerance e and the third column specifies the initial 
constraint tightening Sinit for the stopping condition and the 
relative accuracy requirement for the constraints when using 
optimality conditions. 

Columns four, five and six contain the simulation results. 
The results are obtained by simulating the system with 1000 
randomly chosen initial conditions that are drawn from a 
uniform distribution on X. Column four and five contain the 
mean and max numbers of iterations needed and column six 
presents the average constraint tightening 5 used at termination 
of Algorithm [T] 

We see that the adaptive constraint tightening approach 
gives considerably less iterations for a larger initial tightening. 
However, for more than 10% initial constraint tightening 
((^init = 0.1), the number of iterations is not significantly 
affected. It is remarkable to note that 50% initial constraint 
tightening (Sinit — 0.5) is as efficient as, e.g., 5% ((Jinit = 
0.05) considering that more reductions in the constraint tight- 
ening need to be performed. This indicates early detection 
of infeasibility. We also note that for a suitable choice of 
initial constraint tightening, the average number of iterations 
is reduced significantly. 

VI. Conclusions 

We have equipped the duality-based distributed optimization 
algorithm in |19|, when used in a DMPC context, with a 
stopping condition that guarantees feasibility of the optimiza- 
tion problem and stability and a prespecified performance of 
the closed-loop system. A numerical example is provided that 
shows that the stopping condition can reduce significantly the 
number of iterations needed to achieve these properties. 
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A. Proof for Lemma |2] 

We divide the proof into two parts, the first for x = and 
the second for i 7^ 0. For a; = we have at iteration fc = 
that y*^ = which is the optimal solution. Hence ( |24] | holds 
for k ~Q since all terms are and — A^^_t^ € X. 

Next, we show the result for 5; 7^ 0. Whenever (IT) is feasible 
we have convergence in primal variables ||T9] Theorem 1]. This 
together with the linear relation through which ^ is defined 
( fT4] i gives — ^ z* for T = 0, . . . , — 1 as fc cx). We 
have z* e {1 — S)X and since (1 ~ S)X C A" for every 
S £ (0, 1] this implies that there exists finite fcp such that 
G A" for all k > k^. Equivalent convergence reasoning 
holds for w^. Together this implies that there exists finite fc^ 
such that Pn{x,v'') < 00 and that Pn{x,v'^) — >■ V^{x) for 
all k > kg . Together with convergence in dual function value 
|fT9l Theorem 1] gives that 

D%{x,\\fi'')>PN{x,v'')-et{x) 

holds with finite fc since £* {x) > and e > 0. This concludes 
the proof. □ 

B. Proof for Lemma \3\ 

We introduce y'= = [{i^ {x,6)Y {w'' {x,5)YY , where 
$,^{x,5) and v'^(a;,(5) satisfies the dynamic equations (fT4l l. 
Whenever (|24] | holds we have that ^^'(x, (5) G X and 
v^{Xt5) G U for r — 0,...,A^ — 1. We also introduce 
y* = [(2;*(x,0))^(v*(5,0))^]^. This implies 



i(y'=-y*fH(: 



y'-y*) = 



2(y*)'^Hy* 



< Pn{x, v^-) - V°{x) < D%{x, A^ + et{x) - V^{x) 

< Sin'^fd + etix) 



where the first inequality comes from the first order opti- 
mality condition i27i Theorem 2.2.5] and by definition of 
and Pjv. The second inequality is due to ( l24l i and the 
last inequality follows from Lemma [T] Further, since H = 
blkdiag(Q, . . . ,Q, R, . . . , R) we have for r 0, . . . , iV - 1 
that 

2 



1 






'z*{x,0) 




2 


v';{x,S) 









H 



<^(y^--y*^^ 



H(y^- - y*) 



where H — blkdiag((5, i?), whenever ( l24l i holds. This com- 
pletes the proof. □ 



C. Proof for Lemma |4] 



Since x S but x ^ X'^ we have that V^{x) < 00 and 



V^{x) = 00. Further, from the strong theorem of alternatives 
||28! Section 5.8.2] we know that since V^{x) — 00 for the 
current constraint tightening 6 the dual problem is unbounded. 
Hence there exist \f, fj,j: such that 

6fi^d>D%{x,Xf,,if)-V°ix)>2£tix) (33) 

where Lemma [T] is used in the first inequality. Further, the 
convergence rate in 1,29. Theorem 4.4] for algorithm (IT&i - (ir3b 

is 



Di{x,X*,ti*yD%{x,X\t,'^) < 



2L 



(k + iy 



By inspecting the proof to ||29l Theorem 4.4] (and [29 Lemma 
2.3, Lemma 4.1]) it is concluded that the optimal point 
A* , n* can be changed to any feasible point Xf,fj,j: and the 
convergence result still holds, i.e.. 



D%ix,Xf,tifyD%{x,X\ti'') < 



2L 



(fc + 1) 



That is, there exists a feasible pair (Xf^fij) such that with 
finite fc we have 



D'j,{x, A^ /x'^) > D%ix, Xf,fif) - et{x). 



(34) 



This implies 

mV >^i(a;,A^ yj?(x) 

> D%{x,Xf,^Jif)~V^{x)~ et (x) > et (x) 

where Lemma [T] is used in the first inequality, (|34] | in the sec- 
ond inequality and ( l33b in the final inequality. This completes 
the proof. □ 

D. Proof for Theorem\l\ 

To prove the assertion we need to show that the do loop will 
exit for every x e int(X5^). For every point x G int(X5;) there 
exists 5 e (0, 1) such that G int(X5^). Since int(X5^) C 
X^, we have that V^(Yrj) < °° the optimal solution 
y(_£_,0) satisfies Ay*(^,0) = and Cy*(^,0) < 

d. We create the following vector 



y(x) (!-%*( 



0) 



(35) 



which satisfies 

Ay(x) = Ay*(-^,0)(1-^) =bxi^ =bx (36) 

1 — 1—0 

Cy(x) = Cy*(-^,0)(l-5)<d(l-^). (37) 

1-0 

Hence, by definition ( l23l l of we conclude that for every 
X G int(X5^) there exist 5 G (0,1) such that x G X^. This 
implies that for every x G int(X^) we have that either x G X^ 
for the current constraint tightening 5 € (0,1) or x ^ X^ but 
X G X^^. Thus, from Lemma |2] and Lemma |4] we conclude 
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that either the do loop is terminated or S is reduced and I 
is increased for every x E mt{X'j^) with finite number of 
algorithm iterations k. 

To guarantee that the do loop will terminate for every x e 
int(X5^), we need to show that the conditions in the do loop 
will hold for small enough S and with finite k. That is, we 
need to show that the following two conditions will hold. 

1) For small enough 6, i.e., large enough I, we have that 



< ee*{x) 



(38) 



where 6 = 2 '(5init holds for every algorithm iteration 
k. 

2) For small enough S, i.e., large enough I, the condition 
D%{x,X\^'') > PNiAx + Bv^,^^';)+a£{x,v^) (39) 
with a satisfying ( |25] l holds with finite k whenever 



D%{x,X\fJ-'')>PNix,^^)- 



l + l 



£{x,v^) (40) 



holds. 



We start by showing argument 1 . From the convergence rate 
of the algorithm |fT9l it follows that there exists D > — oo 
such that D^pf{x,\'^,fi'') > D_ for every algorithm iteration 
A: > 0. This is used below where we extend the result from 
ll30l Lemma 1] to handle the presence of equality constraints. 
For algorithm iteration /c > 0, x € int(X5^) and 5 < 6/2 we 
have 

D<D%{x,X\ti') 

= mi iy^Hy + (A'=)^(Ay - hx)+ 
y 2 

+ (MT(Cy - (1 - <5)d) 

< i(y(x))^Hy(x) + (X'^fiAyix) - hx)+ 

+ (/xT(Cy(s) - (1 - <5)d) 

< (1 - + if^' fiCm - (1 - ^)d)+ 

.k\T 



it^ydi6-6) 



<V^ij^^) + {f^'fdi6--6) 



0/ * X l/..fc^Td^ 



1-6' 2^^ 



where the equality is by definition, the second inequality holds 
since any vector y{x) is gives larger value than the infimum, 
the third and fourth inequalities are due to d35l . (l36i and d37] l 
and since {1 — 5) G (0, 1) and the final inequality holds since 
6 < 6/2. This implies that 



which is finite. We denote by U the smallest I such that 6 > 
2~''*5init- Since 6 = 2~'5init this implies that 



6{ti''fd<6 



< 2-'6u 



<2-^+..+i(^^(_^)_^) 



2-'^(5i„it 
(41) 



as ^ — > oo. Especially, with finite I we have that (1381 ) holds 
for every algorithm iteration k. This proves argument 1 . 

Next we prove argument 2. We start by showing for large 
enough but finite I that P'^{Ax + i?i^Ar(x), v^) is finite 
whenever ( |40] | holds. From the definition of P/v and we 
have that {Ax+Bi'n{x),v'^) is finite whenever Pn{x, v^) 
is finite and if A£^'^_^{x,6) S X. For algorithm iteration k 
such that ( l40l i holds we have 



\\A{eN^Ax,6)^z%_,{x,o))r< 



< 



2\\Af 



■\\eN-li=i.S)-Z%^,ix,0)\\H 



(42) 



as ^ — > c» where H = blkdiag(Q, R) and the smallest 
eigenvalue Ainin(-ff) > since H is positive definite. The 
first inequality follows from Cauchy-Schwarz inequality and 
Courant-Fischer-Weyl min-max principle, the second inequal- 
ity comes from Lemma [3] and the third comes from (|411 . 
By definition of we have Az'^_-^^{x,Q) G int(A') which 
through (l42l i implies that A£^'^_^{x,6) G X for some large 
enough by finite I, i.e., small enough 6, and for algorithm 
iteration k such that ( |40l ) holds. 

What is left to show is that ( [39] l holds for every a < 1 — 
2e - K{V2e + V*w)^(\/2e + 1)^ for large enough but finite 
I whenever (l40l i holds. From Lemma |3] and (HTt we know for 
large enough I and any algorithm iteration k such that (l40l) 
holds that 



1 












2 













<6i^ 



k\T, 



H 



l + l 



t{x) 



= 2-'<5i„it(M")^ d + j^t{x) < 2tt{x) 



for any r = 0, . . . , iV — 1, where H = blkdiag(Q, K). Taking 
the square-root and applying the reversed triangle inequality 
gives 



H 



< 



< 2y/ee*{x). 



H 



(43) 
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This implies that 



SAT-l 
<-l 



< 



H 



^JV-1 



H 



H 



'2e] 



so 



H 



^+%/2^)(l + \/2^) 



H 



where we have used ( l43T l, 



ck 
so 



X, 



\J z^Qz + v'^Rv = ^/2£{z, v) and Definition [T] Squaring 
both sides gives through the definition of k that 



< + V2'ef{l + V2'e)H{^^,v^). (44) 

We get for large enough / and for k such that ( |40| | holds that 



i5^(^,A^M'=)> 

>PAr(x,v'=) 



1 



>PAr(i,v'^)-er(x) 
^Pn{Ax + Bv^^,w^^) + {1 
> Pn{Ax + Bv^,v';)+ 

1 - e - + 



2e)^(l 



eK(eo',«o)-^*(^e^-i) 



e{x,v^) 



where the first inequality comes from (|40), the second since 
I > 0, the equality is due to ( |T6] l. the third inequality comes 
from ( |44] |. and the final inequality comes from (l25t . This 
concludes the proof for argument 2. Thus, the do loop will 
terminate with finite I and k. This implies that is defined 
for every x E int(X^), i.e. that dom(i'jv) ^ int(X5^). 
Finally, to show ( |26] l we have that 



V^ix) > D%{x, A^ - 

> Pn{Ax + Bv^,^r^^) - e£*{x) + a£{x,v^) 

> V^iAx + Bv^) + {a- e)£{x, v^) 

where the first inequality comes from Lemma [1] the second 
from (|38] l and (|39] l which obviously hold also for any x £ 
dom(i^Ar), and the third holds since Pn{Ax + Bvq,v'^) > 
Vat (Ax + Bvq) and by definition of £*. This concludes the 
proof. □ 



